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SOLUTIONS  OF  RADIALLY  SYMMETRIC  STRESS 
PROBLEMS  BY  NUMERICAL  METHODS 

Introduction 

Solutions  of  stress  problems  have  for  the  most  part  been 
obtained  with  the  use  of  analytical  methods  which  usually 
involve  series  expansions.  Relatively  little  consideration  has 
been  given  to  the  development  of  the  numerical  methods  that 
are  available  for  the  solution  of  boundary  value  problems. 

The  mathematical  difficulties  in  a  formal  solution  very 
often  make  the  practical  engineering  problem  intractable. 
However  the  numerical  methods,  especially  those  which  in¬ 
volve  displacements,  enable  one  to  obtain  a  solution  for  most 
mixed  boundary  value  problems. 

When  a  body  undergoes  deformation  due  to  some  cause 
stresses  may  be  produced  throughout.  The  nature  of  these 
stresses  is  dependent  on  the  disturbance,  restraints,  geometry 
of  the  body  and  physical  properties  of  the  material.  Utilizing 
the  theory  of  elasticity  for  elastic  problems  and  of  plasticity 
for  inelastic  problems,  a  solution  of  the  stresses  in  a  given 
problem  can  be  obtained  by  solving  under  prescribed  boundary 
conditions  the  differential  equations  which  describe  the  state 
of  deformation  of  the  body.  Since  the  solution  of  a  large  class 
of  practical  problems  is  very  difficult  by  analytical  treatment, 
the  object  of  the  investigation  was  to  develop  satisfactory 
approximate  numerical  methods  for  the  solution  of  stress 
problems.  The  particular  class  of  problems  investigated 
possesses  radial  symmetry. 

The  solution  to  a  stress  problem  is  obtained  at  a  discrete 
number  of  points.  The  system  of  coordinates  used  as  a  frame 
of  reference  and  the  size  of  the  mesh  or  the  network  of  points 
adopted  in  a  solution  depend  on  the  nature  of  the  problem. 
For  many  radially  symmetric  problems  a  cylindrical  co¬ 
ordinate  system  using  a  square  network  of  points  in  an  axial 
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plane  suffices.  However,  for  certain  problems  which  have 
sharp  stress  gradients  due  to  large  curvatures  it  is  expedient 
to  use  an  orthogonal  curvilinear  system  of  coordinates. 

Two  general  methods  were  considered  in  the  development 
of  the  numerical  procedures  for  the  solution  of  stress  prob¬ 
lems.  One  is  the  lattice  analogy  method  in  which  it  is 
imagined  that  the  solid  body  can  be  replaced  by  a  framework 
of  bars.  The  second  scheme  involves  the  well  known  pro¬ 
cedure  of  expressing  a  differential  equation  in  its  finite 
difference  form.  In  both  methods  the  displacements  at  the 
discrete  points  are  considered,  and  a  solution  is  obtained  by 
a  process  of  iteration. 

The  classical  assumptions  made  in  the  theory  of  elasticity 
were  also  adopted  in  the  development  of  the  approximate 
methods  for  the  solution  of  stresses  in  an  elastic  body. 


Development  And  Results  Of  Methods  Used  For 
The  Solution  Op  Elastic  Problems 

The  lattice  analogy  method  of  solving  elastic  stress  prob¬ 
lems  was  considered  first.  The  basic  concepts  in  the 
development  of  the  method  are  simple.  It  is  assumed  that 
the  behaviour  of  the  solid  body  under  load  can  be  approxi¬ 
mated  by  a  framework  of  bars  subjected  to  the  same  load. 
The  pattern  or  arrangement  of  bars  is  arbitrarily  selected, 
but  suitable  to  the  type  of  problem  to  be  solved.  For  the 
chosen  lattice  pattern  the  area  of  each  bar  is  determined 
from  a  consideration  of  the  deformation  of  an  element  of 
the  solid  body  and  the  corresponding  element  of  the  frame¬ 
work,  when  both  are  subjected  to  similar  load  conditions.  Or 
simply,  the  bar  areas  are  so  determined  that  the  deformation 
of  both  elements,  solid  and  lattice,  are  identical  for  certain 
simple  force  systems. 

A  consideration  of  statics  at  a  node,  for  a  system  under 
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arbitrary  displacements,  enables  one  to  develop  the  required 
displacement  equations.  These  equations  for  a  particular 
node  are  related  to  the  force  concentrations  at  the  node  in 
question  and  to  the  displacement  components  at  the  immediate 
surrounding  nodes.  Two  such  equations  exist  at  all  nodal 
points.  The  solution  of  the  set  of  linear  equations  that  occur 
in  a  problem  is  obtained  by  a  process  of  successive  corrections. 

From  the  results  of  the  displacement  components  at  each 
of  the  various  node  points,  forces  in  the  bars  of  the  framework 
are  computed.  By  passing  radial  and  axial  sections  through 
the  center  of  a  node  and  considering  the  bar  forces  on  one 
side  of  a  section,  one  is  able  to  evaluate  an  equivalent  uniform 
stress  distribution  for  the  solid  body.  This  stress  distribution 
is  a  step  function  and  its  approximation  of  the  actual  stress 
distribution  in  the  solid  body  depends  in  part  on  the  coarse¬ 
ness  of  the  net  used  in  solving  the  problem. 

The  lattice  analogy  method  was  found  suitable  for  the 
solution  of  mixed  boundary  value  problems.  However,  the 
method  appears  to  be  limited  due  to  time  consuming  com¬ 
putational  work  to  bodies  whose  geometric  shape  is  formed 
of  straight  radial  and  axial  lines.  Curved  surfaces  in  the 
framework  are  represented  by  a  series  of  steps  which  follow 
in  an  approximate  manner  the  outline  of  the  boundary. 
Problems  solved  with  such  irregular  boundaries  in  the  analogy 
method  showed  high  stress  concentrations  at  the  re-entrant 
corners.  The  solutions  in  these  cases  were  not  representative 
of  the  actual  curved  boundary  problem.  Solutions  with  finer 
subdivisions  in  the  region  of  the  curved  surface  were 
attempted,  but  the  results  did  not  appear  to  justify  the  lengthy 
computational  work  involved. 

In  an  eifort  to  obtain  a  more  suitable  method  of  defining 
a  curved  surface  by  a  discrete  number  of  points  and  of 
sin^lifying  the  solution  of  curved  boundary  problems  to  a 
point  where  a  good  approximation  could  be  obtained  with  a 
coarse  net  in  a  short  period  of  time,  the  lattice  analogy 
method  was  replaced  by  the  method  of  finite  differences. 


4 


The  basic  equations  to  be  integrated  are  the  equations  of 
statics  (1)  which  are  written  in  terms  of  the  radial  and  axial 
component  displacements  u  and  w. 


1  \  1  Kr 

Y2 - Iq  _| - j - ~  o  (a) 

r^/  l—2[x  0r  G 

1  0ev  K, 

V^w  H - 1 - =  0  (^) 
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In  equations  (1),  the  symbol  e^  denotes  the  dilatation,  fx  and  G 
are  the  elastic  constants  and  Kr  and  are  components  of  the 
body  force.  Equations  (1)  are  solved  subject  to  the  pre¬ 
scribed  boundary  conditions  which  may  be  stated  in  terms 
of  stress,  displacements,  or  a  combination  of  both. 

Derivatives  at  interior  node  points  are  defined  by  their 
simplest  finite  difference  expressions.  For  example,  a  first 
derivative  is  specified  by  two  points,  second  derivatives  by 
three,  and  so  forth.  At  the  boundary  fictitious  node^  points 
were  used  and  the  derivatives  were  defined  by  their  simplest 
finite  difference  expressions.  Higher  order  terms  in  the  finite 
difference  expressions  for  a  derivative  were  considered  only 
in  certain  special  cases  of  a  curved  boundary. 

Boundary  conditions  specified  in  terms  of  stress  are 
transformed  into  conditions  of  displacements  by  means  of 
the  stress-strain  and  strain-displacement  relations.  These 
conditions,  when  written  in  finite  difference  form  permit  one 
to  determine  the  displacements  at  the  fictitious  node  points 
in  terms  of  known  displacements  and  surface  stresses. 

Curved  boundary  surfaces  in  the  numerical  solution  of  a 
problem  give  rise  to  difficulties  which  are  inherent  in  the 
adopted  system  of  cylindrical  coordinates.  No  general  method 
appears  to  be  suitable  for  establishing  a  curved  boundary  by 
a  discrete  number  of  points  and  of  defining  the  boundary 
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conditions  thereon.  Actually,  each  problem  has  to  be  con¬ 
sidered  separately  and  the  scheme  used  in  solving  a  problem 
depends  to  a  large  extent  on  the  computer's  judgment  and  ex¬ 
perience.  Four  allied  schemes  for  treating  curved  boundaries 
were  adopted.  In  the  first  scheme,  slopes  at  a  boundary  were 
defined  by  two  points.  In  the  second  scheme,  four  points 
were  used  and  in  the  third,,  six.  The  fourth  scheme  dealt  with 
multiple  displacement  values  at  the  fictitious  boundary  nodes. 

The  various  schemes  were  applied  to  the  same  notch 
problem.  It  was  found  that  all  four  methods  gave  very  nearly 
the  same  values  for  the  axial  stress  at  all  points  of  the  net¬ 
work.  However,  in  the  region  of  the  root  of  the  notch,  neither 
the  results  obtained  from  the  adopted  coarse  or  fine  subdivided 
nets  agreed  with  the  exact  solution  as  given  by  Neuber  in  his 
book  “Theory  of  Notch  Stresses.”  The  disagreement  was 
primarily  due  to  the  coarseness  of  the  net,  which  did  not 
permit  an  accurate  representation  of  the  curved  surface  and 
the  boundary  conditions  thereon.  Further,  the  network  of 
points  in  the  region  of  the  notch  was  not  fine  enough  to 
furnish  an  accurate  description  of  the  stress  variation. 

To  obtain  a  better  solution  for  the  notched  bar  problem 
a  more  convenient  system  of  coordinates  than  the  cylindrical 
system  was  adopted,  namely,  an  orthogonal  curvilinear 
system.  Such  a  system  has  certain  inherent  properties  which 
make  it  desirable  for  the  solution  of  curved  boundary  prob¬ 
lems.  No  difficulties  arise  in  the  delineation  of  a  curved 
surface  by  a  discrete  number  of  points  or  in  the  determination 
of  the  finite  difference  expressions  for  the  boundary  conditions 
in  terms  of  stress  that  apply  to  the  curved  surface. 

The  particular  coordinate  system  used  was  the  ellipsoidal. 
The  network  of  intersecting  traces  of  hyperboloids  and  ellip¬ 
soids  on  an  axial  plane  yielded  a  pattern  which  had  the 
desirable  feature  of  the  finer  and  denser  subdivisions  occuring 
in  the  region  of  high  stress  concentration  at  the  root  of  the 
notch. 
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The  solution  of  the  notched  bar  problem  in  the  curvilinear 
coordinate  system  was  obtained  subject  to  prescribed  bound¬ 
ary  displacements.  The  required  boundary  displacement 
components,  well  away  from  the  root  of  the  notch,  were  taken 
from  a  solution  attained  in  a  cylindrical  coordinate  system. 
Two  solutions  of  the  notched  bar  were  made  with  a  curvi¬ 
linear  pattern  of  points.  One  solution  was  secured  with  a 
network  of  2  units  by  2  units  and  the  second  from  a  network 
of  4  by  4  units.  In  the  latter  case,  the  axial  stresses  along  the 
radial  centerline  were  in  close  agreement  with  the  exact  solu¬ 
tion.  At  the  root  of  the  notch,  this  solution  was  in  error  by 
-20.6  percent;  whereas,  for  the  same  stress,  the  fine  subdivided 
net  in  the  cylindrical  system  was  in  error  by  -39.9  percent. 
By  extrapolating  from  the  two  solutions  made  in  the  curvi¬ 
linear  system  of  coordinates,  it  was  found  that  the  axial  stress 
at  the  root  of  the  notch  differed  by  -9.5  percent  from  the  exact 
solution. 

A  brief  note  regarding  body  forces  in  the  lattice  analogy 
and  finite  difference  methods  is  given  in  the  concluding  part 
of  the  thesis.  The  appendices  treat  the  development  of  the 
lattice  analogy  method  for  the  solution  of  two  and  three- 
dimensional  elastic  problems  with  reference  to  a  rectangular 
system  of  coordinates. 


Development  And  Results  Of  A  Method  Used  For 
The  Solution  Of  Inelastic  Problems 

Although  the  basic  equations  for  a  stress  analysis  in  the 
plastic  range  are  well  known,  mathematical  complexities  have 
limited  the  solutions  of  problems  to  only  a  few  special  cases. 
General  methods  of  integrating  the  fundamental  equations 
in  the  various  theories  of  plasticity  have  not,  as  yet,  been 
developed.  To  obtain  solutions  of  mixed  boundary  value 
problems  a  numerical  method  was  investigated. 

The  basic  stress-strain  equations  used  are  of  the  deforma- 
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tion  type  which  deal  with  small  strains.  These  equations 
are  considered  as  made  up  of  two  parts,  one  part  dealing  with 
the  elastic  deformations  and  the  second  with  the  plastic 
deformations.  The  yield  and  flow  condition  used  is  assumed 
to  be  a  function  only  of  the  second  invariant  of  the  stress 
deviation.  To  define  the  start  of  plastic  yielding  or  to  define 
the  flow  condition  for  perfectly  plastic  materials  the  Huber- 
Mises-Hencky  condition  was  used.  In  the  development,  it  was 
assumed  that  the  material  was  under  a  state  of  gradually 
increasing  load,  and  the  stress  analysis  was  made  at  a 
particular  stage  during  the  loading  process.  Further,  it  was 
also  assumed  that:  (a)  the  axes  of  principal  stress  and  strain 
coincide;  (b)  the  deviations  of  stress  and  strain  are  pro¬ 
portional;  and  (c)  the  intensity  of  stress  is  a  completely 
determined  function  of  the  intensity  of  strain  for  all  complex 
states  of  stress.  The  intensities  of  stress  and  strain,  except 
for  a  constant  factor,  are  equal  to  the  square  root  of  the 
second  invariants  of  the  stress  and  strain  deviations, 
respectively. 

The  basic  differential  equations  to  be  integrated  were 
obtained  by  substituting  stress-strain  relations  of  the  deform¬ 
ation  type  into  the  equations  of  statics.  These  differential 
equations  were  then  written  in  terms  of  displacements. 
Further,  the  equations  were  split  into  two  parts ;  one  part  was 
identical  in  form  to  equations  (1),  as  obtained  in  the  elastic 
theory.  The  second  part  involved  a  function  associated  with 
the  inelastic  behavior. 

The  integration  of  these  differential  equations  subject  to 
prescribed  boundary  conditions  was  obtained  by  first  express¬ 
ing  them  in  their  finite  difference  form  and  then  solving  the 
set  of  linear  algebraic  equations  that  result  in  a  particular 
problem  by  a  process  of  iteration.  During  each  cycle  or  group 
of  cycles  in  the  iteration  process,  the  part  involving  the  in¬ 
elastic  function  was  considered  constant,  and  for  convenience 
it  was  considered  to  be  a  fictitious  body  force  with  components 
in  the  axial  and  radial  directions.  Since  the  fictitious  body 
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forces  and  displacements  are  interrelated,  changes  in  one 
produced  a  change  in  the  other.  Thus,  iterations  were  carried 
out  until  both  the  displacements  and  the  fictitious  body  forces 
became  stabilized. 

The  above  procedure  was  applied  to  a  simple  problem  of 
a  thick  hollow  cylinder  subjected  to  an  internal  pressure  and 
to  two  other  more  complex  problems.  These  were  (a)  a  short 
cylinder  subjected  to  a  tensile  parabolic  stress  distribution 
over  its  ends;  (b)  a  short  cylinder  compressed  between  rough 
rigid  surfaces.  The  solution  for  the  thick  hollow  cylinder 
problem  was  compared  to  existing  solutions  in  which  the 
assumption  of  incompressibility  was  not  made.  A  comparison 
of  the  stresses  was  very  good,  with  variations  of  less  than  1.0 
percent  in  both  the  elastic  and  plastic  ranges. 


Conclusions 

The  lattice  analogy  method  as  developed  is  suitable  for 
the  solution  of  mixed  boundary  value  problems  in  which  the 
geometric  shape  of  the  body  is  formed  of  straight  radial  and 
axial  lines.  Curved  boundaries  impose  a  limitation  on  the 
method,  as  the  framework  can  only  approximate  the  curved 
boundary  by  a  series  of  steps.  At  each  re-entrant  corner  the 
high  stress  concentrations  that  appear  create  a  stress  disturb¬ 
ance  over  a  large  region  of  the  framework.  The  use  of  fine 
subdivided  frameworks  in  the  region  of  the  curved  surface 
to  localize  the  effect  of  each  re-entrant  corner  is  not  desirable 
as  the  computational  work  in  a  solution  becomes  exceedingly 
great  and  the  convergence  very  slow. 

The  finite  difference  method  appears  to  be  suitable  for  the 
solution  of  certain  curved  boundary  problems  when  cylindrical 
coordinates  are  used  as  the  frame  of  reference.  This  system 
of  coordinates  is  adequate  to  determine  a  stress  distribution 
provided  the  geometrical  shape  of  the  solid  body  is  such  that 
no  regions  of  exceedingly  sharp  stress  gradient  occur.  For 
other  problems,  networks  of  fine  subdivisions  in  regions 
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of  large  curvature  are  not  sufficient  to  enable  one  to  obtain 
a  solution  which  would  be  indicative  of  that  of  the  actual 
problem.  In  these  cases  it  is  necessary  to  perform  the 
numerical  work  with  reference  to  a  suitable  orthogonal 
system  of  coordinates. 

The  numerical  method  presented  for  the  solution  of  in¬ 
elastic  problems  deals  with  stress-strain  relations  of  the 
deformation  type.  The  flow  condition  is  considered  to  be  a 
function  only  of  the  second  invariant  of  the  stress  deviation. 
The  basic  equations  to  be  integrated  that  arise  from  these 
fundamental  relations  and  the  equations  of  statics  are  linear 
in  form  when  expressed  in  terms  of  displacement  components. 
This  characteristic  of  linearity  of  the  differential  equations 
was  a  desirable  feature,  as  it  permited  one  to  obtain  a  solution 
of  an  inelastic  problem  by  solving  a  sequence  of  elastic 
problems  through  a  process  of  successive  corrections. 

The  utility  of  the  finite  difference  methods  presented  for 
the  solution  of  certain  elastic  and  inelastic  problems  is  beyond 
question.  However  there  still  remains  the  development  of 
means  of  hastening  the  convergence  of  a  problem  to  a  solution, 
also  of  developing  systematic  methods  of  group  corrections  to 
shorten  the  computational  labour.  Further  investigations 
and  improvements  are  also  required  in  the  definition  of  a 
curved  surface  and  the  boundary  conditions  thereon  by  a 
rectangular  pattern  of  a  discrete  number  of  points. 

The  lattice  analogy  method  in  its  present  form  is  not 
suitable  for  the  solution  of  inelastic  problems.  This  phase 
of  the  method  requires  further  investigation  and  development. 
Also  the  numerical  methods  presented  require  consideration 
with  a  view  to  solving  anisotropic  problems  for  both  the 
elastic  and  inelastic  behavior. 


